## This file uses the combined dataset for analysis to run all regressions ##
## Note that you will need to have obtained the Leip datasets with county-level turnout, ##
## otherwise your dataset will be missing the dependent variable ("turn") ##
## Created by Meredith Dost and last run 8/19/2025 ##

# load packages
library(plm)
library(lmtest)
library(clusterSEs)
library(fixest)
library(tidyverse)
library(broom)
library("DIDmultiplegtDYN") 

# set working directory
# setwd("")

# load in data--all U.S. counties
burd <- read.csv("final_dataset_for_analysis_replication.csv")
## all border counties w/in 100 miles of border ##
burd100 <- subset(burd, MINDIST<=100)


##### MODELS WITH 3 INDICES #####
turn.bothidx <- as.formula(turn ~ regidx + turnidx + negsum + exp | year + fips)
turn.bothidx_full <- as.formula(turn ~ regidx + turnidx + negsum + exp
                                + higheligib + exp:higheligib +
                                  + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod <- feols(turn.bothidx, burd, cluster = "fips")
turn.bothidx_full.mod <- feols(turn.bothidx_full, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100 <- feols(turn.bothidx, burd100, cluster = "fips")
turn.bothidx_full.mod.100 <- feols(turn.bothidx_full, burd100, cluster = "fips")

##### MODELS WITH 3 INDICES + MEDICAID POPULATION #####
turn.bothidx.mp <- as.formula(turn ~ regidx + turnidx + negsum*pcteligib + exp  | year + fips)
turn.bothidx_full.mp <- as.formula(turn ~ regidx + turnidx + negsum*pcteligib+exp 
                                   + exp:higheligib +
                                     + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod.mp <- feols(turn.bothidx.mp, burd, cluster = "fips")
turn.bothidx_full.mod.mp <- feols(turn.bothidx_full.mp, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.mp <- feols(turn.bothidx.mp, burd100, cluster = "fips")
turn.bothidx_full.mod.100.mp <- feols(turn.bothidx_full.mp, burd100, cluster = "fips")

##### CTY POP WEIGHTED MODELS WITH 3 INDICES #####
### running models ###
## among all US counties
turn.bothidx.mod.wt <- feols(turn.bothidx, burd, cluster = "fips", weights=~lpop)
turn.bothidx_full.mod.wt <- feols(turn.bothidx_full, burd, cluster = "fips", weights=~lpop)

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.wt <- feols(turn.bothidx, burd100, cluster = "fips", weights=~lpop)
turn.bothidx_full.mod.100.wt <- feols(turn.bothidx_full, burd100, cluster = "fips", weights=~lpop)


##### MODELS WITHOUT ELECTION BURDEN #####
turn.med <- as.formula(turn ~ negsum + exp | year + fips)
turn.med_full <- as.formula(turn ~ negsum + exp
                            + higheligib + exp:higheligib +
                              + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.med.mod <- feols(turn.med, burd, cluster = "fips")
turn.med_full.mod <- feols(turn.med_full, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.med.mod.100 <- feols(turn.med, burd100, cluster = "fips")
turn.med_full.mod.100 <- feols(turn.med_full, burd100, cluster = "fips")

##### MODELS WITH COMBINED ELECTION BURDEN INDEX #####
turn.bothidx.e <- as.formula(turn ~ cov + negsum + exp | year + fips)
turn.bothidx_full.e <- as.formula(turn ~ cov + negsum + exp
                                  + higheligib + exp:higheligib +
                                    + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod.e <- feols(turn.bothidx.e, burd, cluster = "fips")
turn.bothidx_full.mod.e <- feols(turn.bothidx_full.e, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.e <- feols(turn.bothidx.e, burd100, cluster = "fips")
turn.bothidx_full.mod.100.e <- feols(turn.bothidx_full.e, burd100, cluster = "fips")


##### MODELS WITH COVI MEASURE #####
# subset dataset to presidential years and remove DC to be consistent with COVI analyses
burd.p <- subset(burd, year %in% c(2012,2016,2020))
burd100.p <- subset(burd100, year %in% c(2012,2016,2020))
burd.p <- subset(burd.p, state!="DC")
burd100.p <- subset(burd100.p, state!="DC")

# get sample sizes for state and county
length(unique(burd.p$state))
length(unique(burd.p$fips))
length(unique(burd100.p$state))
length(unique(burd100.p$fips))

turn.bothidx.cv <- as.formula(turn ~ FinalCOVI + negsum + exp | year + fips)
turn.bothidx_full.cv <- as.formula(turn ~ FinalCOVI + negsum + exp
                                   + higheligib + exp:higheligib +
                                     + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod.cv <- feols(turn.bothidx.cv, burd.p, cluster = "fips")
turn.bothidx_full.mod.cv <- feols(turn.bothidx_full.cv, burd.p, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.cv <- feols(turn.bothidx.cv, burd100.p, cluster = "fips")
turn.bothidx_full.mod.100.cv <- feols(turn.bothidx_full.cv, burd100.p, cluster = "fips")


##### MODELS -- PRESIDENTIAL YEARS #####
### running models ###
## among all US counties
turn.bothidx_full.mod.p <- feols(turn.bothidx_full, burd.p, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx_full.mod.100.p <- feols(turn.bothidx_full, burd100.p, cluster = "fips")

##### MODELS -- MIDTERM YEARS #####
# subset dataset to presidential years
burd.m <- subset(burd, year %in% c(2010,2014,2018))
burd100.m <- subset(burd100, year %in% c(2010,2014,2018))
### running models ###
## among all US counties
turn.bothidx_full.mod.m <- feols(turn.bothidx_full, burd.m, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx_full.mod.100.m <- feols(turn.bothidx_full, burd100.m, cluster = "fips")


##### MODELS WITH FOX INDEX #####
turn.bothidx.f <- as.formula(turn ~ regidx + turnidx + foxneg + exp | year + fips)
turn.bothidx_full.f <- as.formula(turn ~ regidx + turnidx + foxneg + exp
                                  + higheligib + exp:higheligib +
                                    + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod.f <- feols(turn.bothidx.f, burd, cluster = "fips")
turn.bothidx_full.mod.f <- feols(turn.bothidx_full.f, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.f <- feols(turn.bothidx.f, burd100, cluster = "fips")
turn.bothidx_full.mod.100.f <- feols(turn.bothidx_full.f, burd100, cluster = "fips")

##### MODELS WITH CONDENSED INDEX AB3 #####
turn.bothidx.r3 <- as.formula(turn ~ regidx + turnidx + AB3 + exp | year + fips)
turn.bothidx_full.r3 <- as.formula(turn ~ regidx + turnidx + AB3 + exp
                                   + higheligib + exp:higheligib +
                                     + pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + swing + sen*gub  | year + fips)

### running models ###
## among all US counties
turn.bothidx.mod.r3 <- feols(turn.bothidx.r3, burd, cluster = "fips")
turn.bothidx_full.mod.r3 <- feols(turn.bothidx_full.r3, burd, cluster = "fips")

## among US counties w/in 100 miles of border
turn.bothidx.mod.100.r3 <- feols(turn.bothidx.r3, burd100, cluster = "fips")
turn.bothidx_full.mod.100.r3 <- feols(turn.bothidx_full.r3, burd100, cluster = "fips")

## saving out
save.image("mods_output.RData")


#### MODELS WITH ESTIMATORS ROBUST TO HETEROGENEOUS TREATMENT EFFECTS #####
## uses package DIDmultiplegtDYN version 1.0.15
burd100$sen_gub <- burd100$sen*burd100$gub
burd100$exp_higheligib <- burd100$exp*burd100$higheligib

covars.use <- c("regidx","turnidx","exp",'higheligib',"pctwhite","pcths","pct65plus","l_medincome",'l_cvap','demshare','swing',"sen_gub","exp_higheligib")

mod.tot <- did_multiplegt_dyn(df = burd100, outcome = "turn", group = "fips", time = "year", treatment = "negsum",
                              continuous=1, bootstrap = 500, controls=c(covars.use), effects=1)

mod.tot.nocovar <- did_multiplegt_dyn(df = burd100, outcome = "turn", group = "fips", time = "year", treatment = "negsum",
                                      continuous=1,bootstrap=500,controls=c("regidx","turnidx","exp"))

#cleaning up R environment
rm(list=ls()[!ls() %in% c("mod.tot","mod.tot.nocovar","burd100")])

## saving out
save.image("bootstraps_mods_output.RData")

#cleaning up R environment
rm(list=ls()[!ls() %in% c("burd100")])

####### HETEROGENEOUS TREATMENT TESTING ####### 
# Build models for the residualized outcomes
out_resid <- lm(turn ~ regidx + turnidx + exp*higheligib + 
                  +pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + 
                  swing + sen * gub  + factor(year) + factor(fips), data = burd100)

trt_resid <- lm(negsum ~ regidx + turnidx  + exp*higheligib + 
                  +pctwhite + pcths + pct65plus + l_medincome + l_cvap + demshare + 
                  swing + sen * gub  + factor(year) + factor(fips), data = burd100)

fpe_weights <- burd100 %>%
  mutate(treatment_resid = residuals(trt_resid)) %>%
  mutate(treatment_weight = treatment_resid / sum(treatment_resid^2)) %>%
  mutate(mean_treat = mean(negsum),
         diff = negsum - mean_treat)

fpe_weights %>%
  summarize(beta = sum(turn * (treatment_resid / sum(treatment_resid^2))))

fpe_weights %>%
  summarize(twfe_beta = sum(turn * treatment_weight))

fpe_weights %>%
  summarize(total_weights = sum(treatment_weight))

# Add residuals to data with weights
fpe_weights <- fpe_weights %>%
  mutate(out_resid = residuals(out_resid))

# run models to evaluate slope
check_slopes <- lm(out_resid ~ treatment_resid * negsum,
                   data = fpe_weights)

summary(check_slopes)


## saving out
save.image("heterogeneouschecks.RData")
#cleaning up R environment
rm(list=ls()[!ls() %in% c("")])

####### CCES 2016 DATA regressions #######
# set working directory
# setwd("")

# read in data
d2 <- read.csv("other_data/cces_burden_final.csv")
d.all <- d2
# subset to adults < 65 years old
d2 <- subset(d.all, age<=64)

#### cces among ALL adults < 65 ####
out_d2 <- feglm(
  vv ~  negsum1612*highmed+exp_alt*higheligib + turnidx + regidx +
    blacknh+othernh+hisp+ faminc + age + educ + female | region,
  data = d2, weights=~commonweight_vv_post, family="binomial", cluster = "fips_state")

#### cces among adults < 65 who are NOT Medicaid recipients ####
out_not <- feglm(
  vv ~  negsum1612*highmed+exp_alt*higheligib + turnidx + regidx +
    blacknh+othernh+hisp+ faminc + age + educ + female | region,
  data = subset(d2, med==0), weights=~commonweight_vv_post, family="binomial", cluster = "fips_state")

#### cces among adults < 65 who *ARE* Medicaid recipients ####
out <- feglm(
  vv ~  negsum1612*highmed+ exp_alt*higheligib + turnidx + regidx +
    blacknh+othernh+hisp+ faminc + age + educ + female |region ,
  data = subset(d2,med==1), weights=~commonweight_vv_post, family="binomial", cluster = "fips_state")

#### cces among ALL-AGED adults ####
out_d.all <- feglm(
  vv ~  negsum1612*highmed+exp_alt*higheligib + turnidx + regidx +
    blacknh+othernh+hisp+ faminc + age + educ + female | region,
  data = d.all, weights=~commonweight_vv_post, family="binomial", cluster = "fips_state")


## saving out
#save.image("cces_output.RData")
